clearvars; clc; close all;
%%
dt_pre = load('pre_regulation.mat');
dt_reg = load('regulation.mat');

dz = dt_pre.dz;
z = dt_pre.z;
beta = dt_pre.beta;
R_pre = dt_pre.R;
R_reg = dt_reg.R;
z_min = dt_pre.z_min;

Q_pre = exp(dt_pre.qstar);
Q_reg = exp(dt_reg.qstar);
%%

%q1_start = exp(z_min + beta*R_pre-1);
amin = 3;
amax = 23;
N = 100;
x1 = linspace(amin,amax,N)';
dx1 = x1(2) - x1(1);
dt_reg.f = dt_reg.f/sum(dt_reg.f(Q_reg>=amin & Q_reg<=amax)*dz);
f_reg = dt_reg.f*dz/dx1;
[bch1,dist_reg1] = gp_dist(Q_reg,f_reg,x1);

err =  normrnd(0,0.04,[5e4,1]);
q_err = exp(dt_reg.qstar + err);
[~,dist_err] = gp_dist(q_err,f_reg,x1);

data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab.csv');
data_orig.z         = data_orig.assets./data_orig.brnumber_hh; % Create additional variables

% Starting with correct values: -------------------------------------------
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;
fpost = ones(length(apost),1);
drl = 1/length(apost);
f_rl = fpost*drl/dx1;
[~,dist_rl] = gp_dist(apost,f_rl,x1);

figure
hold on
plot(bch1,dist_rl,'+','LineWidth',1)
% plot(bch1,dist_reg1,'--','LineWidth',1)
plot(bch1,dist_err,'LineWidth',1)
xlabel('Asset Size')
xlim([amin amax])
 set(gca, 'yScale', 'log')
  set(gca, 'xScale', 'log')
y = ylim; % current y-axis limits
plot(exp(dt_reg.zl(1) + beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
plot(exp(dt_reg.zh(1) + beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
legend('Data','Model')
hold off



%%

%q1_start = exp(z_min + beta*R_pre-1);
amin = 30;
amax = 70;
N = 40;
x1 = linspace(amin,amax,N)';
dx1 = x1(2) - x1(1);
dt_reg.f = dt_reg.f/sum(dt_reg.f(Q_reg>=amin & Q_reg<=amax)*dz);
f_reg = dt_reg.f*dz/dx1;
[bch1,dist_reg1] = gp_dist(Q_reg,f_reg,x1);

err =  normrnd(0,0.04,[5e4,1]);
q_err = exp(dt_reg.qstar + err);
[~,dist_err] = gp_dist(q_err,f_reg,x1);

data_orig           = readtable('bhc_and_commercial_hh_min10_branches_for_matlab.csv');
data_orig.z         = data_orig.assets./data_orig.brnumber_hh; % Create additional variables

% Starting with correct values: -------------------------------------------
cond        = data_orig.assets >= amin & data_orig.assets <= amax & (data_orig.post == 0 | data_orig.post ==2); %& data_orig.datequarter_num <= 2017.5;
data        = data_orig(cond,:);
apre  = data(data.post == 0,:).assets;
apost = data(data.post == 2,:).assets;
fpost = ones(length(apost),1);
drl = 1/length(apost);
f_rl = fpost*drl/dx1;
[~,dist_rl] = gp_dist(apost,f_rl,x1);

figure
hold on
plot(bch1,dist_rl,'+','LineWidth',1)
% plot(bch1,dist_reg1,'--','LineWidth',1)
plot(bch1,dist_err,'LineWidth',1)
xlabel('Asset Size')
xlim([amin amax])
 set(gca, 'yScale', 'log')
  set(gca, 'xScale', 'log')
y = ylim; % current y-axis limits
plot(exp(dt_reg.zl(2) + beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
plot(exp(dt_reg.zh(2) + beta*dt_reg.R - 1)*[1 1], [y(1) y(2)],'k-','LineWidth',.5);
legend('Data','Model')
hold off